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Abstract. Quantum cascade lasers can be modeled within a hierarchy of differ- 
ent approaches: Standard rate equations for the electron densities in the levels, 
semiclassical Boltzmann equation for the microscopic distribution functions, and 
quantum kinetics including the coherent evolution between the states. Ifere we 
present a quantum transport approach based on nonequilibrium Green functions. 
This allows for quantitative simulations of the transport and optical gain of the 
device. The division of the current density in two terms shows that semiclassical 
transitions are likely to dominate the transport for the prototype device of Sirtori 
et al. but not for a recent THz-laser with only a few layers per period. The many 
particle effects are extremely dependent on the design of the heterostructure, and 
for the case considered here, inclusion of electron-electron interaction at the Ifartree 
Fock level, provides a sizable change in absorption but imparts only a minor shift 
of the gain peak. 

1 Introduction 

Since the first realization of a quantum cascade laser (QCL) in 1994 [1] these 
semiconductor heterostructures have become important devices in the in- 
frared regime operating up to room temperature [2] . Lasing in the THz-region 
was also achieved recently [3], opening a new window for applications. The 
standard devices contain an injector region guiding the electrons to the upper 
laser level in an active region where the optical transitions occur between a 
few discrete levels. A frequently studied prototype is the sample in [4]. Differ- 
ent designs are interminiband-QCLs [5,6], as well as QCLs without injector 
regions [7] or containing only four barriers per period like the staircase-laser 
[8] and a recent THz-QCL [9]. 

The modeling of quantum cascade lasers was first performed on the basis 
of rate equations [10] for the electron dynamics in the active region. It was 
assumed that the electrons reach the upper laser level with the rate J/e, 
where J is the current density and e < is the electron charge. A necessary 
condition for inversion is that the scattering rate I/t^^i from the upper to 
the lower laser level is smaller than the out scattering 1/r/ from the lower 
laser level. Optimizing these scattering rates by a sophisticated choice of well 
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and barrier widths in the active region, QCLs with high performance could 

be designed. While typically scattering with optical phonons is considered 
to be the main mechanism for the scattering rates [11,12], electron-electron 
scattering has also been treated [13,14]. The influence of a magnetic field has 
been studied by these rate equations in [15]. While these rate equations for 
the electron densities rii [in units l/cm^j for the levels i average over the 
momentum k in the in- plane direction, the distribution functions /i(k) can 
be taken into account employing Monte-Carlo (MC) simulations [16-18]. 

If one includes the injector region in the simulation and imposes periodic 
boundary conditions (a good approximation as typical devices have approx- 
imately 30 periods each containing an active region and an injector region) 
a full simulation of QCL-devices can be performed. Such an approach was 
performed almost simultaneously on the basis of rate equations [19], MC- 
simulations [20] and a quantum transport model [21] obtaining good results 
for the current-voltage characteristic of a prototype device [4] . 

In this article we want to show, in how far quantum effects affect the 
transport and gain behavior and address the question if simple semiclassical 
models such as rate equations or MC-simulations are applicable. In particular 
we demonstrate (i) that the current can be calculated in a quantum transport 
model, (ii) how this relates to semiclassical approaches, and (iii) discuss the 
implications of many particle corrections on the gain spectra. 

2 Current in Quantum Transport 

In order to describe the quantum cascade laser we start by defining a set 
of single particle basis states ^a{z)e'^''^ / VA. Here k, r are two dimensional 
vectors in the x, y plane perpendicular to the growth direction z and A is the 
normalization area. The functions ^a{z) refiect the layer sequence of the QCL 
structure and may be chosen as cnergy-eigenstates or Wannier states (see the 
discussion in [22]). Then the Hamilton operator reads in second quantization: 



where all terms connecting different k-indices (i.e. breaking the translational 
invariance of the structures) have been included in fl^scatt- The spin index s 
yields an additional factor 2 for the current and the gain as we assume that 
all states are spin degenerate and no spin transitions occur. 

Much information is contained in the (reduced) density matrix 




(1) 



a,/5,k,s 



/5ak,/3k' 




(2) 



in particular the occupation probabilities are given by the diagonal elements 
/a(k) = Pakak- The average current density (in the ^-direction) is evaluated 
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by the temporal evolution of the position operator z 
^ ' ^ , ' ^ , ' 

— J^ — </scatt 

where V denotes the normalization volume. Let us first consider the current 
Jq. For an arbitrary choice of the basis we may write 

Jo = 2(for Spin):^ ^ ^ W0,c«(k)pc<k,/3k (4) 

c«/3k 

where 

7 

is an anti-hermitian matrix. If the wave functions {z) are chosen real, which 
is typical for bound states, W^j,a(k) becomes real and Jq is determined by 
the non-diagonal elements of Pak,/3k- 
For a scattering part of the form 

^scatt = X] Oak,7k'(t)a]jk(0a7k'(0 1 (6) 
a7k,k',s 

which contains only pairs of electronic annihilation and creation operators, 
we obtain 



'Oak,0k'(i)^/37 ~ ^a/30/3k,7k' tt^k' ) • (7) 



Jscatt - ^ ^ \^«k 



ak 7/3k' 

In the case of phonon scattering Oak. 7k' contains phonon annihilation and 
creation operators and thus phonon-assisted density matrices [23] determine 
Jscatt for this scattering process. 

To evaluate the density matrices we perform the perturbation expansion 
within the formalism of nonequilibrium Green functions [24-26] similar to 
[27] . The key quantities are the lesser and retarded Green function 

G'a„a.(k;ii,i2) = i(aj,^k(*2)aaik(ti)) (8) 

GLta.(k;ii,t2) = -i6>(ti -t2)(a„,k(ii)4^k(*2) +al,k(*2)aaik(ii)) (9) 

where the time dependence is taken in the Heisenberg picture. The lesser 
Green function refers to the electron density and it becomes the density 
matrix Paik,a2k(0 = ^J]&.\t^i)l\ for equal times. In the stationary state 
considered here the Green functions only depend on the time difference t = 
t\ — t-2 and we introduce the energy E as the Fourier conjugate of t: 

G{t2 + 1, h) = I ^G(£)e-'^*/'^ (10) 
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This provides us with the Dyson equation 

{E - H°(k) - E'^'^iyL, E)) G''°'(k, E)^\ (11) 
and the Kcldysh relation 

G<(k,£;) = G'''^*(k,£)i:<(k,£;)G''^^(k,E) (12) 

where capital bold symbols represent matrices in a/?. Together with the func- 
tionals i^{G} for the self-energies this provides a self consistent set of equa- 
tions which can be solved numerically. Although the Green functions are 
diagonal in k, the expression (7) for Jscatt can be evaluated by Eq. (21) as 
derived in the appendix. 

Here we use self-energies in self-consistent Born approximation for im- 
purity, interface roughness, and phonon scattering, applying the following 
approximations: (i) The k-dependence of the scattering matrix elements is 
neglected, (ii) It is assumed that 17 is diagonal and depends only on the diag- 
onal elements of G in the basis of Wannier functions. The scattering matrix 
elements are evaluated for a typical momentum transfer assuming an inter- 
face roughness with average height of 0.28/\/27r nm ^ and a correlation length 
of 10 nm. The impurity scattering was estimated by an effective scattering 
rate 7imp/?i- Electron-electron interaction is included within the mean field 
approximation. See [22] for further details. 

3 The Current- Voltage Chciracteristic 

We perform our calculation using a basis of Wannier functions. These func- 
tions are shown in Fig. 1 for zero bias and an operating field of 220 mV per 
period for the sample used in [4] . While the spatial structure of the Wannier 
functions does not change with bias, their energetic position is affected both 
by the external field and the mean field which is evaluated self consistently. 
From Fig. 1 we see, that the mean field almost vanishes at operating con- 
ditions as the electrons are mainly in the injector region where the doping 
is also located. The energy levels of the Wannier functions bunch at the op- 
erating field indicating the strong coupling between the functions enabling 
transport through the structure. 

In Fig. 2 the current-voltage characteristic is shown for different doping 
densities A^u per period. The theoretical result exhibits a monotonic increase 
of the current density with doping, showing that the mean field has no dra- 
matic influence on the transport behavior in these structures. We find good 
agreement with the experimental data except for Nu = 3.5 x 10^^/cm-^ and 
Nd = 5.3 X 10^^/cm^ where the experiment exhibits a significantly higher 
bias drop. The difference is of the same order as the bias of a test structure 

^ In the calculations performed in [22] a factor 2-w was lacking in the program, 
which can be compensated by the reduction in the roughness height. 
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z [nm] z [iimj 

Fig. 1. Conduction band ofTset including mean field potential and Wannier func- 
tions for two different electric fields for the sample of [4] 




Fig. 2. Current- voltage characteristics for different doping densities for the struc- 
tures of [28] at r = 77 K. (a) Experimental data from M. Giehler (PDI Berlin) 
where the thin line refers to a structure grown without the QCL structure, thus 
providing an estimate for contact effects, (b) Theoretical result for 7imp = 5 meV 
using the bias U — NFd, where the QCL structure consists of TV = 30 periods 



containing only contact layers and the waveguides, albeit it is not clear why 
this additional bias drop is only present in some samples. 
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Fig. 3. Contributions to the current density for tlie samples of [4] (a) and [9] (b). 
7imp = was used in the calculations 



4 Comparing Jq and Jgcatt 

In Fig. 3(a,b) we show the different contributions to the current evaluated 
for the structures of [4] and [9], respectively. Both current-field relations are 
in reasonable agreement with the respective experimental results. (The data 
of [9] only extends to Fd w 70 mV, therefore there is no verification of 
the current peak.) While Jscatt dominates the behavior in Fig. 3a, both the 
contributions of Jq and Jscatt are important in Fig. 3b. 

In the following we want to discuss the role of the two current contribu- 
tions with respect to the use of semiclassical approaches: As the expressions 
for Jo and Jscatt are invariant to unitary transformations of the basis states, 
they can be evaluated in arbitrary basis sets. A special basis set is given 
by the energy eigenstates f^iz) obtained by diagonalizing H° (including the 
mean field), which will be used in the following argumentation. 

The semiclassical theories used in [19,20] imply that the density matrix is 
diagonal in the energy eigenstates, i.e., /O^k/j'k — '^^^'/^(k). In this basis the 
diagonal elements of W^ij,' in Eq. (5) vanish and thus Jq becomes zero in the 
semiclassical approximation. 

In the semiclassical approximation the Green functions in the basis of 
energy eigenstates are given by 

G-y-<i^(k, E) « TniS^^,6iE - E^) 
G<,(k, E) « 2^i(5^,,5(£; - £;^)/^(k) 

Then we find from Eq. (21) 

-^scatt ~ "jT^ X! { 2^il^^' ' 



(13) 
(14) 



i;^) + i/4k)i:-*^(k,i;^) 



-i:<^(k,i;^) + i/^(k)i:-*(k,E^) 



(15) 
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In semiclassical approximation the self-energies are related to the scattering 
probabilities -R^'k'->/:ik as follows 

and the quantities Z'^^ and iJ''***^ contain an additional factor z^'n'. This 
provides us with 

•^scatt ~ -j^ -R/ik^/i'k'(^A('A»' ~ ^Mm) (^^) 

which is the semiclassical expression for the current density^. Therefore the 
entire current is contained in Jgcatt in the semiclassical approximation. 

For the special case of the structure considered in [4], the density matrix 
is approximately diagonal in the basis of energy eigenstates^ implying that 
Jo — > and Jscatt is well approximated by the semiclassical expression (16). 
This expectation is supported by Fig. 3a and the findings of [20]. In contrast, 
Jo is an important contribution for the THz-laser of [9], which contains only 
4 barriers per period, see Fig. 3b. Therefore it is questionable if semiclassical 
approaches work here. 



5 Gain and Absorption Spectra 

The general evaluation of gain spectra within the quantum transport model 
used here was described in detail in [29] . The key idea is to evaluate the com- 
plex susceptibility xi^) which is related to the optical absorption coefficient 
at a frequency uj via [30] 

a(u) = , 

c ub 

where ns is the background refractive index and c is the speed of light. 
Figure 4 shows the gain spectrum for the sample of [4]. At zero current we 
find strong absorption due to transitions in the active region, which vanishes 
already for small currents as the carriers are transfered to the injector region. 
Pronounced gain around Tio; = 130 meV sets in for current densities of several 
kA/cm^. The height and width of the gain spectrum is in good agreement 
with the findings of [31]. 

^ Terms of the form z^ti/0^k_p'k'G^'^'(k)0^/k/ have been neglected for 7^ ^ 

here. Their implication is not clear yet. 
^ Note that Jo also vanishes if the density matrix is diagonal in a basis of real 

states, which are not energy eigenstates. But then Jscatt no longer corresponds 

to the semiclassical result. 



(17) 
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Fig. 4. Gain spectrum for tlic 
14 16 sample of [4] using 7imp = 0. 

(Prom [29]) 



Within the semiclassical approximation the susceptibihty is given by [32] 

where is the FWHM of 3{G[^°)(k, E)} and p^^ = ez^^ is the dipolc matrix 
element. In [22] it was shown that this semiclassical approach gives reasonable 
results compared with the quantum model for the structure of [4]. 

In these approaches the influence of electron-electron interaction was to- 
tally neglected. Here we study the influence of many-particle corrections 
within the Hartree Fock approximation on the gain spectrum. The suscepti- 
bility is decomposed by 

X{uj) = 2(for spin) ^ p^^X!^,/:i(k, w) 

where the susceptibility functions Xi/,/i(k, w) between the eigenstates f and 
H are determined by the equation 

p.^(k) iuik) - Uk)) = 

n{u- e,(k) + e^(k) + i(r^ + r,)/2) x.^(k, w) 
+ iUik) - Uik)) 2 V ^J" ^) E^-M(k',^) (19) 





(/.(k)-^(k)) Y,x.^{^',t)v 



v V /i /i 

k k' 



Equation (19) reduces, in the equilibrium case, with only two isolated sub- 
bands of idealized quantum well subbands where phenomenological dephasing 
characterizes the broadening, to Eq. (5) of [33]. The bare Coulomb interaction 
and renormalized energies which appear above are given by 

^ ( k - k' ) - J d.d.>;(.)^.(.) ,,^,^^n,_,,/| C(^O^M^O 
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Fig. 5. Absorption spectrum 
for the sample from Fig. 2 
with Nd = 3.8 X 10" /cm^ at 
an operating field of Fd = 0.2 
V. The dashed line gives the 
result from Eq. (18). The full 
line includes Coulomb correc- 
tions according to Eq. (19) 



E/.(k')^(7/k')- 

The second term on the right-hand side of (19) gives rise to the depolar- 
ization shift [34,35] , while the last term (exchange contribution) is analogous 
to the excitonic coupling term in interband transitions [36]. Figure 5 shows 
the absorption spectra. The inclusion of many-particle corrections yields a 
blue shift of about 5 meV for the low frequency absorption peak and a slight 
red shift for the gain peak around 130 meV. 

6 Discussion 

The impact of quantum effects on transport and gain in quantum cascade 
lasers have been examined. In the evaluation of the current two different 
terms, Jq oc [Ho,z] and Jscatt oc [-ffscatt,^] appear. In the semiclassical ap- 
proximation, where the density matrix is assumed to be diagonal in the basis 
of energy eigenstates, Jscatt carries the entire current. Our quantum trans- 
port calculations show that Jscatt dominates the behavior for the prototype 
sample of [4], which has been frequently studied, thus justifying the semiclas- 
sical approaches in [19,20]. On the other hand the current Jo, resulting from 
nondiagonal elements in the density matrix, shows strong contributions for 
the THz-laser of [9]. Nevertheless, it is not clear by now in how far the as- 
sumption of diagonal self-energies in the Wannier basis affects this behavior. 
Ongoing work is focused towards the inclusion of the full matrix structure in 
the self- energies. 

The many particle effects are extremely dependent on the structure, since 
its design determines the actual electronic overlap and subband occupation. 
For the case considered here, the gain spectra are hardly modified by the 
electron-electron interactions within the Hartree-Fock approximations, while 
a significant depolarization shift occurs for the low frequency absorption. 
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Appendix 

Similarly to [22] the scattering current (7) can be evaluated in the following 
way: We define the contour-ordered Green function (superscript c) 

-P'ak(n,T2) = -iTc{(Oak,/3k'(ri)Z/3T,a^k'(n)a),k(T2) 
7/3k' 

- 2;a/30/3k,7k'(Tl)a^k'(n)a]jk('^2))} • 

They are evaluated in the Dirac representation (with index D) 



f'ak(n,T2) = ^ -it,(e/<i-**-"(-){6fk,/3k'(n)^/37<k'(n)af^(T2) 

7/3k' 

-^«/30^k,7k'(n)a:°,,(n)af^(r2)}; 
The lowest order non- vanishing terms of the expansion gives 

7/3k' Se 

(0„k,/3k'(n)^/37G7/(k'; Ti, r)65k',ek(T)G^°„(k; T, T2) 

- ^a/30/3k,7k'(Tl)G;°5(k'; Ti, r)65k',ek(T)GJ°Jk; T, T2) 



with the bare Green functions G^°^(k;ri,r) = —iTc{{a^]^{Ti)d^^{T))}. In 
order to be consistent with the perturbation expansion in the Green functions, 
further terms arc taken into account, which replace the bare Green functions 
by the full Green functions. Then we find 



f'ak(n,T2) « ^E/'^^ [r-(k;n,r)G^,Jk;T,T2) 



X) ^a/3^^e(k; Ti, T)G^_„(k; T, Ta) 
/3 



with 



i:^,^(k;Ti,T) = (Oak,/3k'(Tl)2;37G^,5(k';Tl,r)05k',.k(T)) (20) 

7/3ik' 

where the averaging refers to the phonon bath for phouon scattering. Thus, 
in Born approximation the self energies are given the the usual function- 
als for the self-energies -^(G) where the Green functions G are replaced by 
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Z • G in matrix notation. Using Langreth rules and changing to the energy 
representation t) can be inserted in Eq. (7) yielding the final expression 

■^catt ^ -^ak(^' 

ak 

_2e ^ f dE f 

ak (21) 

e 

- J2 [^|e(k, E)Gt^{K E) + Sf:{^^, ^)G< Jk, E)] } 
to evaluate Jscatt- 
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